Kinetic energy density functionals from the Airy gas, with an apphcation to 
the atomization kinetic energies of molecules 



Lucian A. Constantin and Adricnn Ruzsinszky 
Department of Physics and Quantum Theory Group, Tulane University, New Orleans, LA 70118 

(Dated: June 27, 2009) 

We construct and study several semilocal density functional approximations for the positive Kohn- 
Sham kinetic energy density. These functionals fit the kinetic energy density of the Airy gas and 
they can be accurate for integrated kinetic energies of atoms, molecules, jellium clusters and jellium 
surfaces. We find that these functionals are the most accurate ones for atomization kinetic energies 
of molecules and for fragmentation of jellium clusters. We also report that local and semilocal 
kinetic energy functionals can show "binding" when the density of a spin unrestricted Kohn-Sham 
calculation is used. 
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I. INTRODUCTION 



The positive Kohn-Sham (KS)^ kinetic energy (KE) 
density of noninteracting electrons 



1 ^ 



(1) 



is an exact functional of the occupied orbitals 
Density functional approximations to the noninteract- 
ing kinetic energy r5[n|,nj = J drT(r) can simplify and 
speed up by orders of magnitude any KS self-consistent 
calculation^. (Here n|(r) and ni{r) are the spin densi- 
ties.) However, in spite of important and hard work done 
in this direction^,, no actual approximation has reached 
chemical accuracy. 

The simplest model of an edge electron gas is the Airy 
gas, where any electron feels a linear effective potential^, 
and thus the normalized one-particle eigenfunctions are 
proportional to the Airy function. The effective finite- 
linear-potential model gives remarkably good results for 
the jellium surface problem^. However, the KE den- 
sity derived in this approximation^ does not recover the 
correct second-order gradient expansion KE density^i^ 
and has an unphysical oscillating behavior in the limit 
of slow density variations'^, being a poor approximation 
for atoms'l. 

The positive KE density of the Airy gas was studied 
by Vitos et. ali^, and they derived a generalized gra- 
dient approximation (GGA) density functional for T(r). 
(This approximation is denoted in this paper by VJKS 
GGA.) They showed that the poor behavior of the ki- 
netic energy density derived in the linearized-potential 
approximation^ is mainly due to a Laplacian term that 
arises naturally in the Airy gas model. Thus, the Lapla- 
cian term, even if it integrates to zero and does not affect 
the integrated KE, is an important tool in developing 
density functionals not only for the KE but also for the 
exchange-correlation (xc) energ y '^i'^ . The VJKS GGA 
KE density functional fits the Airy gas KE density and 
is a good model for the KE density of the jellium sur- 
faces, but for atoms and molecules it diverges to —oo at 



the nuclei, due to the behavior of the Laplacian term. 
The integrated kinetic energies are at a Thomas-Fermiii 
level of accuracy, reducing considerably the error of the 
linearized-potential approximation^. 

A jellium surface is the simplest model of a metal- 
lic surface. Self-consistent local-spin-density (LSD) 
calculations'^ for this model provided early evidence that 
density functionals may work. But wavefunction-based 
methods, like Fermi hypernetted chain'- and Diffusion 
Monte Carlo (DMC) of Reff^i predicted low-density sur- 
face xc energies about 40% larger than those from LSD. 
Recent refined DMC estimatesi^i, and calculations in the 
random phase approximatio n'^i^" and beyond i t^'i^^i^'^ , 
agree with the popular xc semilocal density functionals, 
showing that the jellium surface can not only be accu- 
rately described in the context of density functional the- 
ory, but can also be an important model used to develop 
new density functionals. 

The exchange energy density of the Airy gas^i^liS^ 
and the xc jellium surface energie o^^i^^ were employed 
in the construction of accurate xc GGA's for solids. (See 
Refsj^i^^iS^). A simple xc GGA functional depends only 
on spin densities and their gradients and can not de- 
scribe accurately both solids and atoms^. However, a 
Laplacian-level xc meta-GGAH, that depends nontriv- 
ially on spin densities and their gradients and Laplacians, 
can be accurate for atoms, molecules, solids and surfaces. 

In this paper, we derive several GGA KE function- 
als from the Airy gas and jellium surfaces and we find 
them accurate for atomization KE energies of molecules 
and for fragmentation of jellium clusters. Our function- 
als, constructed similarly to that of RefJ^, recover the 
second-order gradient expansion of the integrated KE, 
have the right behavior of the KE density in the tail of 
the density, and fit the kinetic energy density of the Airy 
gas. 

The paper is organized as follows. In section [Hi we 
construct our KE functionals. In section Hill we test the 
functionals for atoms, jellium clusters, jellium surfaces 
and molecules. In section ITVl we summarize our conclu- 
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II. LAPLACIAN-DEPENDENT GGA KINETIC 
ENERGY FUNCTIONALS 

The positive kinetic energy density of the local Airy 
gas (LAG) iiM 



LAG 



{z)^~^n{z)v,}f{z) + ^W^n{z), (2) 



where Veff{z) is the effective potential and n{z) is the 
density of the Airy gas. (Unless otherwise stated, atomic 
units are used throughout, i.e., = fi = — 1.) Al- 
ternatively, Eq. ^ can be writteni^ using the Thomas- 
Fermi kinetic energy density r^^ — (3/10)(37r^)^/'^n^/'^: 



LAG 



(^) 



TF 



{z)P{z) + -V2n(z), 



where 



P(z) 



2Bz 

(3712)2/3^(^)2/3' 



(4) 



and B is the slope of the linear effective potential. P{z) 
is a smooth function of the reduced density gradient 



3(r) = |Vn(r)|/[2fcF(r)n(r)], 



(5) 



where kp{r) = (37r2n(r))^/'' is the Fermi wavevector. 
(The dimensionless density gradient s(r) measures the 
variation of the density over a Fermi wavelength Xp = 
2TT/kp.) Thus, Vitos et. al^ proposed the following 
GGA KE density functional 



^VJKS 



(r) 



TF/ 



(r)P^-'^^(s(r)) + -V2n(r), (6) 



where 



P 



V.IKS 



1 + 0.8944s2 0.0431s^ 
l + 0.6511s2-f 0.0431s4 



(7) 



fits P{z) for the Airy gas model. Eq. ([6]) recovers 
the exact KE density of the von Weizsacker functional^ 
Vn|2/(8n) — (5/3)r"^^s2 for an exponentially decaying 
density (see KeiM-), but for a slowly- varying density be- 
haves as r'^^(l + 0.2433s2 + 0{s^)) + ^V^n{z) and vi- 
olates the second-order gradient expansion (GE2) of the 
KE density^i^ 



GE2 



TF 



(1 



(8) 



Let us consider the following arbitrary partition of Eq. 
([3]) for the Airy gas model 



LAG 



(^) 



TF 



{z)F{z,P) + pV^n. 



Eqs. dS]) and dH) give 

Fiz,P) = P{z) + 



(1/5) - /3]V2n(z) 



rTF 



(9) 
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F(z, /?) is a smooth function of the reduced gradient s 
for any /3 > 1/8, and it can be accurately approximated 
by the following expression 



1 (ai + 5/27)s2 



1 + ais2 + ass^ -I- 



40/3-5 



(11) 

where 01,02,03,04 and 05 are parameters that depend 
on (3. Eq. pT|) recovers the terms 1 -I- (5/27)s2 for 
a slowly-varying density, but the second-order gradient 
expansion of the KE density additionally requires that 
(3 = 1/6. In the tail, where the density decays exponen- 
tially, Eqs. (Ill) and (fTT|) give the correct KE density of 
the von Weizsacker functional. 

When /3 = 1/5, F(z, /3 = 1/5) = P{z) and we define a 



(3) GGA (Ai) similar with the one in Refif 



TF I 



{r)F^''{s{r),f3 = 1/5) + -V'n{r), (12) 



where the fitting parameters are shown in Table HI 



TABLE I: Parameters of the enhancement factor F^^{s,f3) 
for various GGAs. 





Ai-GGA 


Ai-GGA 


A0.185-GGA 


ai 


1.122609 


1.301786 


1.293576 


12 


0.900085 


3.715282 


2.161116 


13 


-0.227373 


0.343244 


-0.144896 


04 


0.014177 


0.032663 


0.025505 


0.5 


0.731298 


2.393929 


1.444659 



When /3 = 1/6, we define a GGA (A^) that recovers 
the second-order gradient expansion KE density 



T^^(r) 



TF/ 



(r)P^^^(s(r),/3 = l/6) + -VMr), (13) 



where the fitting parameters are shown in Table [H 

The Airy gas is the simplest edge electron gas and does 
not include curvature corrections that are present at the 
edge surfaces (see Fig. 2 of Ref4) . Thus in order to find 
an optimum value of P for jellium surfaces, let us define 
the quality factor (similarly to Refs.^2, and^S) 



m 



dr 



approx 



{r,l3)-T{r)\/ J drrir), (14) 



where t°-pp'^°^ is an approximation of the positive Kohn- 
Sham KE density r. [See Eq. (P)]. We apply the 
quality factor to jellium surfaces using numerical LSD 
Kohn-Sham orbitals and densitiesi^i^. The integration 
was done from Zmin = — 2.75Af to Zmax = 2Ai?, where 
Ai? — 2tt /kp is the bulk Fermi wavelength , for several 
values of bulk parameter Vg. (Here = {%tt/AY/^ /kp is 
the radius of a sphere which contains on average one elec- 
tron, and kp is the bulk Fermi wavevector.) For r°'PP'^°^ 
we use Eqs. (O and (fTT|) . Thus, for values of /3 between 
0.15 and 0.22, we accurately fit F{z,(3) with the Fade 
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0.15 0.16 0.17 0.18 0.19 0.2 0.21 0.22 
P 

FIG. 1: The quality factor i5(/3) versus /?, for r""'"'"'' given 
by Eqs. (O and for the jellium surfaces with bulk pa- 

rameters Ts = 2, 3, and 4. We use LSD KS orbitals and 
densitiesi^i^. 
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FIG. 2: The exact function F{z,/3) shown with points 
{F{z, 13), s{z)) for some discrete z, and parametrized function 
F'^'^is, P) shown with lines for /3 = 1/5, 1/6 and 0.185, versus 
the reduced gradient s, for the Airy gas model. Also shown 



are the enhancement factor (1 -|- 5/27s ) and P 



'is) 



approximation of Eq. pT|) . and we calculate S{P). Fig. 
[T] shows that 5{P) is minimum for (3 « 0.185 for semi- 
infinite jellium surfaces with = 2, 3, and 4. 

So from our jellium surface analysis we define the fol- 
lowing GGA (AO. 185) that also fits the kinetic energy 
density of the Airy gas 

^Ao.i85(j.) ^ T^F(r)F'^'^{s{r),l3 = 0.185)-f 0.185V2n(r), 

(15) 

where the fitting parameters are shown in Table |T1 

In Fig. [2] we show the exact function F{z,P) and the 
fitting function F'^^{s,f3) versus the scaled density gra- 
dient s, for /3 = 1/5,1/6 and 0.185 respectively. Up to 
s = 3, the exact functions F and the parametrized ones 
can not be distinguished. (We note that s values big- 
ger than 3 are found in the tail of an atom or molecule, 
where the electron density is negligible.) P^'^^^{s) over- 
estimates P{z) = F{z, (3 = 1/5) until s « 3 and underes- 
timates P{z) for 3 < s < 10. 

Far from the edge of the Airy gas, the density has 
Friedel oscillation^. These oscillations are well described 
by the kinetic energy density of the linear potential 
approximation^ that in the slowly- varying density regime 
reduces toiS 

H„ TF ^ 5 (Vn)2 1 (Vn)2 2(37r^)i^ 

r =T H sm ; \ . 

72 n 12 n ^ |Vn| ' 

(16) 

The third term represents quantum oscillations and has 
an unphysical behavior when Vn — > 0. In Fig. [3] we show 
T — versus Ci for a slowly- varying Airy gas density. 
The edge is at C = 0. [C = (2i?)^/^z is the scaled spatial 
coordinate for the Airy gas.] The Friedel oscillations are 
well described by Eq. ([16]). But even if r"^^ — t'^^ is the 
worst kinetic energy density shown in the figure, its inte- 
gration over a period of the Friedel oscillations is almost 
exact. Thus t^'s , that behaves as t^^'^ in this limit, is 



0.003 




-0.0005 I ' ' ' ' 1 

-21 -20.8 -20.6 -20.4 -20.2 -20 



FIG. 3: r — r versus for the Airy gas model. The edge 
is at = 0. The integrations of t — r^^ over the complete 
Friedel oscillation shown in figure are: T^'''^' - TJ^ = 8.30 x 

10-^ Tt^ - Tj^ = 9.58 X 10-*, r/^ - rj^ = 8.26 x 10"*, 

yAO.185 _ j,TF ^ g gg ^ j,Hn _ j,TF ^ g gg ^ 

VJKS GGA, not shown in the figure, gives an integrated value 
of 10.09 X 10"". 



the best approximation for the integrated KE, whereas 
^Un gjygg worst integrated KE. 



III. TESTS OF OUR GGA KINETIC ENERGY 
FUNCTIONALS 

In this section we test our functionals for various 
systems. In the calculations we use the spin-scaling 
relation^ 

r<,(M,r) = (l/2)T([n = 2n,],r), (17) 
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where is the density of the electrons with spin cr. 
{<J =T or I) 



A. Integrated kinetic energies of atoms, jellium 
clusters and jellium surfaces 

In Table M we show the accuracy of Tj^, T/'-'^-^, 
rj.GE2^ rj^GEi^ rj^Ak^ j,^^ ^^oms, jel- 

lium clusters and jellium surfaces (similarly as Table I 
of Ref>^). The error displayed in this table is 

Error = ^ "m.a.r.e.atoms" + | "m. a. r.e. clusters" 

+i "m.a.r.e.LDM(N = 8)" , (18) 

where "m.a.r.e. atoms" is the mean absolute relative 
error (m.a.r.e.) of the integrated kinetic energy of 50 
atoms and ions (listed in Ref.— )^ "m.a.r.e. clusters" is 
the m.a.r.e. of 2e~, 8e~, 18e~, 20e", 34e~, 40e~, 58e", 
92e~, and 106e~ neutral spherical jellium clusters (with 
bulk parameter = 3.93 which corresponds to Na), and 
"m.a.r.e. LDM(N=8)" is the m.a.r.e. of the KE of N=8 
jellium spheres for rg — 2, 4, and 6, calculated in the 
liquid drop model^l (LDM) 

j,LDM ^ (3/io)fc2^iV + a^N^^Hnrl (19) 

where kp is the bulk Fermi wavevector, and CTs is the 
surface KE. The exact LDM value is computed with the 
exact as (using LSD orbitals). Because the relative er- 
rors of surface kinetic energies are much larger than those 
of the atoms and spherical jellium clusters, we use the 
LDM approach for calculating the jellium surface KE er- 
rors (as in Refii^); LDM gives m.a.r.e. comparable to 
that of atoms and clusters (see Table |TT| . We use ana- 
lytic Hartree-Fock densities and orbitals'^^ for atoms and 
ions, and numerical Kohn-Sham densities and orbitals for 
jellium clusters (using the optimized potential method 
(OPM)34) jeUium semi-infinite surfaces (using LSD 
xc potential). 

^VJKS ^ ^A^ ^ ^A^ g^j-^^ ^Ao.i85 g.^.^ constructed to model 
the KE density of the Airy gas, but only t^'s recovers the 
second-order gradient expansion of the KE density. The 
difference between t^'s and t^^^^ is given mainly by 
the quality of fitting the function P{z) of Eq. (Sec 
Fig. m) includes effects of density variations 

near jellium surfaces because of our optimization of the 

Laplacian coefficient. In Table [II] we see that Ts is 
very accurate (comparable with the fourth-order gradient 
expansion) for jellium systems and gives an overall error 

smaller than T^^ and TY''^^. Tf°-^^^ is accurate for 
atoms and gives an overall error comparable with the 
fourth-order gradient expansion one ( see also Table 1 of 
ReLi^). 



TABLE II: Mean absolute relative error (m.a.r.e.) of kinetic 
energies of 50 atoms and ions (see Ref .— ) , of neutral spherical 
jellium Na clusters (2e~, 8e~, 18e~, 20e~, 34e~, 40e~, 58e~, 
92e~, and 106e~) and of jellium surfaces (with = 2, = 4, 
and r-g = 6) incorporated into the liquid drop model (LDM) 
for a jellium sphere with N=8 electrons (see Eq. (|19fl '). Also 
shown is the total error given by Eq. psp . 





m.a.r.e. 

atoms 


m.a.r.e. 

clusters 


m.a.r.e. 

LDM(N=8) 


Error (Eq. ITSll) 




0.0842 


0.0439 


0.0810 


0.0733 


rpVJKS 


0.0399 


0.0465 


0.0754 


0.0504 


rpGE2 


0.0112 


0.0099 


0.0330 


0.016 


rpGE4 


0.0251 


0.0176 


0.0170 


0.0212 


rpH 

^ s 


0.0626 


0.0566 


0.0879 


0.067 


^ s 


0.0789 


0.0154 


0.0177 


0.048 


71^0.185 

s 


0.0083 


0.0249 


0.0535 


0.024 



B. Integrated atomization kinetic energy for a set 
of molecules 

In Table Unl we present the atomization kinetic energies 

A- 

for the molecules used in Refs.— i^. We observe that Ts 
keeps the right sign for all the molecules and has practi- 
cally the same mean absolute error as the Thomas-Fermi 
functional. In Ref^S^ it was shown that the Thomas- 
Fermi KE functional gives better atomization kinetic en- 
ergies than all the other tested semilocal functionals. 
jiAo.iss jg accurate for atoms and molecules, and gives 
the smallest mean absolute error for the atomization en- 
ergies presented in Table IIIII We also show that the 
PBE-like semilocal functional of Rcf.'^'*, whose parame- 
ters are fitted to atoms, works worse than the Thomas- 
Fermi functional and all the semilocal functionals derived 
from the Airy gas. 

C. Binding energy of the N2 molecule 

In Fig. d] we show the binding energy of the N2 
molecule as a function of the distance between the nu- 
clei. We use a spin unrestricted Hartree-Fock calculation 
in which the spin symmetry breaks close to the Hartree- 
Fock equilibrium bond length. This helps the functionals 
to show an equilibrium length close to the exact. Figure 
|4] is in accord with the values for the iV2 molecule listed 
in Table mH all the semilocal functionals presented in the 
figure give bigger atomization kinetic energies than the 
exact calculation, thus showing a minimum in the total 
energy calculated with the Hartree-Fock density. 

The unrestricted solution becomes energetically lower 
beyond the Coulson-Fisher point^S than the energy of 
the restricted solution, and spin symmetry breaking for 
the N2 molecule can be achieved by mixing the high- 
est occupied and lowest unoccupied orbitals**^. For a 
spin-restricted calculation, the orbital-free KE function- 
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TABLE III: Integrated atomization kinetic energy ( KE 
atoms - KE molecule, in a. u.) for the set of molecules used 
in Refs.— The kinetic energies were calculated using the 
PROAIMV code with Kohn-Sham orbitals given by the Gaus- 
sian 2000 code (with the uncontracted 6 — 311-|-G(3d/, 2p) ba- 
sis set, Becke 1988 exchange functional^^, and Perdew-Wang 
correlation functional^). The last line shows the mean abso- 
lute errors (m.a.e.). Here TJ^ is the the GGA of Ref.— with 
the parameters k — 0.8438 and /i — 0.2319. 





'T-'exact 

s 




rpVJKS 

^ s 


rpGE2 


a 1 

^ s 


^ S 


rpAO.185 


rpTW 


Ha 


-0.150 


-0.097 


-0.086 


-0.114 


-0.080 


-0.114 


-0.096 


-0.108 


HF 


-0.185 


-0.305 


-0.369 


-0.186 


-0.422 


-0.173 


-0.311 


-0.226 


HaO 


-0.304 


-0.308 


-0.455 


-0.136 


-0.531 


-0.169 


-0.369 


-0.209 


CH4 


-0.601 


-0.737 


-0.907 


-0.571 


-0.972 


-0.618 


-0.813 


-0.649 


NH3 


-0.397 


-0.231 


-0.457 


-0.060 


-0.525 


-0.165 


-0.364 


-0.155 


CO 


-0.298 


-0.323 


-0.580 


-0.085 


-0.678 


-0.181 


-0.456 


-0.203 


F2 


-0.053 


0.128 


0.013 


0.282 


-0.050 


0.269 


0.093 


0.223 


HCN 


-0.340 


-0.1835 


-0.539 


0.079 


-0.644 


-0.097 


-0.399 


-0.071 


N2 


-0.158 


0.344 


-0.046 


0.565 


-0.134 


0.321 


0.069 


0.412 


CN 


-0.431 


-0.215 


-0.539 


0.005 


-0.631 


-0.168 


-0.424 


-0.129 


NO 


-0.268 


0.092 


-0.215 


0.330 


-0.313 


0.176 


-0.094 


0.198 


O2 


-0.100 


0.106 


-0.089 


0.335 


-0.177 


0.286 


0.030 


0.239 


m.a.e. 




0.177 


0.133 


0.311 


0.172 


0.224 


0.116 


0.232 




-6e-05 I ' ' ' ' ^ 
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z(a.u.) 



FIG. 5: t(z) - r"'""'°"(2), where r"'""'°" is r^i r^^ , 
^Ao.i85^ and r*^^^ respectively, versus z, for a jellium sur- 
face of bulk parameter Ts = 3. The surface is at 2 = 0, 
the jellium is at z < and the vacuum is at 2: > 0. 
The surface kinetic energies are: crj^'"^' = —703 erg/cm^, 

at^ = -869 erg/cm^ crf^ = -690 erg/cm^, af°-^^^ = 
-788 erg/cm^ and crf^^ = -762 erg/cm^ VJKS GGA, 
not plotted in the figure, gives aY^'^ ~ —837 erg/cm^. 
(Ihartree/bohr^ = 1.557 x lO^'erg/cm^) We use LSD KS 
orbitals and densitie a^^'^^ . 



CtJ 0.2 



TF ^'aO\\ 


/A0.185 exact ^. 




■' meta-GGA - 



0.6 0.8 1 1.2 1.4 1.6 1.8 
R(N-N) [A*'] 

FIG. 4: Binding energy ( AE = E molecule - E atoms, in 
a.u.) as a function of N-N distance for the N2 molecule using 
a nonrestricted Hartree-Fock calculation (with uncontracted 
6-311-|-G(3df,2p) basis set). The curve "meta-GGA" is the 
binding energy given by the Laplacian-level KE meta-GGA 
of Ref.—. The Hartree-Fock density was used as input for 
orbital free KE functionals. lA = 1. 8897a. w.. 



als listed in Table II do not show an equilibrium point, 
thus the spin-breaking symmetrj«^i^'^ and the spin- 
scaling relations^^ play an important role in describing 
stretched molecules, and they need to be taken into ac- 
count in the orbital-free codes. 



D. Tests of the kinetic energy density 

In Fig. [5] we show the kinetic energy density of our 
functionals at a jellium surface. Though t^^.iss -j^g^ 
smallest overall error, r^e gives the most accurate sur- 
face kinetic energy because it is accurate near the surface 
and it can almost exactly damp the Friedel oscillations 
far from the surface (see Fig. [3]). 

In Fig. [6] we show the kinetic energy densities of our 
functionals for the 2e~ Na jellium cluster. Here the exact 
curve is the von Weizsacker^* KE density. We see that 
all three functionals (r^^s , t"*6, and r'^"-^^^) recover the 
exact curve in the tail of the density, as expected. 



E. Large-Z asymptotic behavior 

The non-interacting kinetic energy of the neutral 
atoms has the following asymptotic expansion^i^: 

= coZ^/^ + ciZ2 +c2Z5/^ (20) 

where Z is the atomic number, and cq = 0.768745, 
ci = -1/2, and C2 = 0.2699. In Ref4^ the authors pro- 
pose an accurate method to extract these coefficients for 
any KE functional. In Table IIVI we present the large-Z 
asymptotic behavior of our functionals. All the function- 
als listed in Table IIVI are exact for systems with uniform 
density, such that we expect that they have the exact 
Thomas-Fermi coefficient cq = 0.768745. (Similarly with 
Ref4^, we do not have enough data points to extract 
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FIG. 6: Kinetic energy density versus radial distance r, for 
the 2e~ jellium cluster (with bulk parameter rs = 3.93). The 
area under the curve is the kinetic energy: yj^^"'^' = 0.114 a.u., 



T 

^ s 



: 0.098 a.u., Ts 



: 0.121 a.u. andT,^ 



VJKS GGA, not shown in the figure, gives 
a.u.. 



= 0.108 a.u.. 
■•^■^^ = 0.101 



Co accurately.) T^^o iss and T/-'^-^, the functionals that 
give the most accurate atomization kinetic energies, have 
reasonable large-Z asymptotic behaviors. 



TABLE IV: The coefficients of the asymptotic expansion of 
Eq. (|20|l for several semilocal functionals. The fitting method 
is the same as in Ref.iS,. We use OPM^ densities. 





CO 


Cl 


C2 


Exact 


0.768745 


-0.500000 


0.269900 


rpGE2 


0.768745 


-0.536197 


0.335992 




0.768745 


-0.507979 


0.291815 




0.768745 


-0.532065 


0.229370 


rp 6 
^ S 


0.768745 


-0.439745 


0.392152 


rpAO.185 


0.768745 


-0.491080 


0.302999 


rpVJKS 

^ s 


0.768745 


-0.507589 


0.225358 



FIG. 7: Error of disintegration KE (DKE'^'^'''^' - DKE^pp''"'') 
for 273 configurations described by Eq. (|21|l . The first point 
iV = 1, corresponds to (106e") (92e") -I- (8e") -1-3 x (2e"), 
and the last point N=273 corresponds to (106e~) ^ 53 x 
(2e~). We use OPM-KS orbitals and densities. Mean ab- 
solute errors are: m.a.e'^'^* = 0.247, m.a.e'*'^ = 0.218, 
m.a.e'^'^^ = 0.128, m.a.e"^^ = 0.092, Tn.a..e'^°-^^^ = 0.039, and 
m.a.e^e = 0.196. VJKS GGA, not shown in the figure, has 
m.a.e^'^^ = 0.165. 



In Fig. El we show DKE°'''''=* - dKE^pp™'' for 273 pro- 
cesses described by Eq. (|2T|) . for several KE functionals. 



We see that our functional Ty^o iss is very accurate, im- 
proving over TJ^ for all the configurations. T^*^ is close 
to, but better than the fourth-order gradient expansion 
T^^*. Overall, this figure agrees well with the atomiza- 
tion KE of molecules reported in Table IIIIl showing an 
important link between jellium spheres and molecules. 
These results and the liquid drop model (see Eq. (17) of 
Ref.— ) suggest that the TF functional gives a good bal- 
ance between jellium surface KE and jellium curvature 
KE. This balance, that is important in atomization and 
disintegration processes, is improved by the A0.185-GGA 
functional. 



Fragmentation of jellium clusters 



IV. CONCLUSIONS 



Let us consider the disintegration of the 106e~ neu- 
tral spherical jellium Na cluster into smaller closed-shell 
jellium spheres: 



(106e") — > ni(92e") + n2(58e") + n3(40e") 
+n5(20e") + n(i{18e^) + n7(8e") 



n4(34e" 



In this paper we have studied several semilocal KE den- 
sity functionals derived from the Airy gas. These func- 
tionals, that depend trivially on the Laplacian of the den- 
sity, do not satisfy several important constraints. Their 
kinetic energy densities are not always positive, and they 
implicitly violate the important constraint t"'PP'^°'^ > 



"■8(2e J^l) (}iere is the von Weizsacker KE density) and diverge 



where ni, ...,ns are positive integers, and 92rii + 58n2 + 
4071-3 -I- 34^4 -1- 20n5 -I- ISng -f Snj + 2ns = 106. We define 
the disintegration KE as 



DKE = KE of initial cluster 



KE of the fragments. 

(22) 



to —oo at the nucleus of an atom. 

However, such functionals can be accurate for the in- 
tegrated KE of jellium surfaces and jellium clusters (e.g. 

T^^), and of atoms and molecules (e.g. 2^"-^*^), when 
we use realistic densities ( from KS calculations). More 
importantly, they are the most accurate KE density func- 
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tionals, to our knowledge, for the integrated atomization 
kinetic energies of molecules and for the fragmentation of 
jellium clusters. These functional may also be useful for 
quasi-realistic densities (e.g. a superposition of free-atom 
Kohn-Sham densities), but they are not accurate enough 
for orbital- free calculations. 

We have also presented a spin-unrestricted Hartree- 
Fock calculation for the stretched A'2 molecule that ex- 
plains the N2 atomization kinetic energies displayed in 



Table IIIIl and that shows equilibrium lengths for many 
semilocal functionals. Thus, this work suggests that the 
spin-symmetry breaking and the spin scaling relations 
can be important tools in orbital- free approaches. 
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